R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

source("tianfengRwrappers.R")

heatmapTFs <- topTFs %>% group_by(seuratCluster) %>% slice_max(n = 2000, order_by = RelativeActivity)

top_regulonActivity <- regulonActivity[sapply(heatmapTFs$Regulon, function(e) {which(rownames(regulonActivity) == e)}), ]
# annotation_row <-  data.frame(cluster = factor(top5$seuratCluster), row.names = make.names(top5$Regulon,TRUE))
genes_to_show <- c(as.character(top5$Regulon), "PRDM6_extended (53g)","PRDM6 (51g)")


##获得基因和细胞聚类信息
cluster_info <- colnames(regulonActivity)

##筛选矩阵为要画热图的基因
mat <- top_regulonActivity

#获得要展示的基因在热图中的位置信息
gene_pos <- match(genes_to_show, rownames(mat))
row_anno <- rowAnnotation(gene=anno_mark(at=gene_pos,labels = genes_to_show))

col <- colors_list[c(4,2,3,1,5)]
names(col) <- cluster_info
top_anno <- HeatmapAnnotation(cluster=anno_block(gp=gpar(fill=col),labels = cluster_info,
                                                 labels_gp = gpar(cex=1,col='black'))) ## 顶端的cluster注释


col_fun <-  colorRamp2(c(-2, 0, 2), c("#1E90FF", "white", "#ff2121")) #颜色

svg("ds2TFs.svg",height = 6,width = 10)
Heatmap(mat, cluster_rows = FALSE, cluster_columns = FALSE, 
        show_column_names = FALSE, show_row_names = FALSE,
        column_split = cluster_info, top_annotation = top_anno,  
        column_title = NULL, right_annotation = row_anno, 
        heatmap_legend_param = list(
          title='Regulon Activity', title_position='leftcenter-rot'), col = col_fun)

Warning: The input is a data frame, convert it to the matrix.

dev.off()

null device 
          1 

ds1

regulonActivity <- read.csv("./ds1_SCENIC/regulonActivity.csv",row.names = 1)
topTFs <- read.csv("./ds1_SCENIC/topRegulators.csv")

regulonActivity <- regulonActivity[,levels(topTFs$seuratCluster)] #调换列顺序
# colnames(regulonActivity)[colnames(regulonActivity)=="T.cell"] <- "T cell"

top5 <- topTFs %>% group_by(seuratCluster) %>% slice_max(n = 5, order_by = RelativeActivity)

#根据已知的top5更新行顺序
top_regulonActivity <- regulonActivity[sapply(top5$Regulon, function(e) {which(rownames(regulonActivity) == e)}), ]


annotation_row <-  data.frame(cluster = factor(top5$seuratCluster), row.names = make.names(top5$Regulon,TRUE)) #make.names用来生成不冲突的行

pheatmap(top_regulonActivity, breaks = unique(c(seq(-2,2, length=400))), 
                 color = colorRampPalette(c("#1E90FF", "white", "#ff2121"))(400),
                border_color = NA, cluster_rows = FALSE, cluster_cols = FALSE,
                main = "regulonActivity",angle_col = 45, show_rownames = T)

# ggsave("TFs_activity.png",device = png,width = 10,height = 8,plot = TFs_heatmap)

heatmapTFs <- topTFs %>% group_by(seuratCluster) %>% slice_max(n = 2000, order_by = RelativeActivity)

top_regulonActivity <- regulonActivity[sapply(heatmapTFs$Regulon, function(e) {which(rownames(regulonActivity) == e)}), ]
# annotation_row <-  data.frame(cluster = factor(top5$seuratCluster), row.names = make.names(top5$Regulon,TRUE))
genes_to_show <- c(as.character(top5$Regulon),"DLX6_extended (43g)","DLX2 (18g)")##对齐
##获得基因和细胞聚类信息
cluster_info <- colnames(regulonActivity)

##筛选矩阵为要画热图的基因
mat <- top_regulonActivity

#获得要展示的基因在热图中的位置信息
gene_pos <- match(genes_to_show, rownames(mat))
row_anno <- rowAnnotation(gene=anno_mark(at=gene_pos,labels = genes_to_show))

col <- colors_list[c(2,1,5,8)]
names(col) <- cluster_info
top_anno <- HeatmapAnnotation(cluster=anno_block(gp=gpar(fill=col),labels = cluster_info,
                                                 labels_gp = gpar(cex=1,col='black'))) ## 顶端的cluster注释


col_fun <-  colorRamp2(c(-2, 0, 2), c("#1E90FF", "white", "#ff2121")) #颜色

svg("ds1TFs.svg",height = 6,width = 10)
Heatmap(mat, cluster_rows = FALSE, cluster_columns = FALSE, 
        show_column_names = FALSE, show_row_names = FALSE,
        column_split = cluster_info, top_annotation = top_anno,  
        column_title = NULL, right_annotation = row_anno, 
        heatmap_legend_param = list(
          title='Regulon Activity', title_position='leftcenter-rot'), col = col_fun)

Warning: The input is a data frame, convert it to the matrix.

dev.off()

null device 
          1 

library(Seurat)
library(SCENIC)
library(AUCell)
library(RcisTarget)
library(SCopeLoomR)
library(dplyr)
library(foreach)

lamb <- function(s1,s2)
{
  paste0(s2,s1)
}

fileloc <- scenicOptions@fileNames

temp <- lapply(fileloc[["output"]],lamb,"ds2_SCENIC/") %>% as.character() %>% as.matrix()
rownames(temp) <- rownames(fileloc[["output"]])
colnames(temp) <- "fileName"
fileloc[["output"]] <- temp

temp <- lapply(fileloc[["int"]],lamb,"ds2_SCENIC/") %>% as.character() %>% as.matrix()
rownames(temp) <- rownames(fileloc[["int"]])
colnames(temp) <- "fileName"
fileloc[["int"]] <- temp

scenicOptions@fileNames <- fileloc
scenicOptions@settings[["tSNE_filePrefix"]] <- "ds2_SCENIC/int/tSNE"

regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
regulonActivity_byCellType <- sapply(
  split(rownames(cellInfo), cellInfo$CellType),
  function(cells) rowMeans(getAUC(regulonAUC)[, cells])
)

regulonActivity_byCellType_Scaled <- t(scale(t(regulonActivity_byCellType), center = T, scale = T))
topRegulators <- reshape2::melt(regulonActivity_byCellType_Scaled)
colnames(topRegulators) <- c("Regulon", "seuratCluster", "RelativeActivity")
topRegulators <- topRegulators[which(topRegulators$RelativeActivity > 0), ]

aucellApp <- plotTsne_AUCellApp(scenicOptions, logMat)

savedSelections <- shiny::runApp(aucellApp)

Listening on http://127.0.0.1:7492
DLX5_extended (13g) threshold replaced by 0.15
DLX6_extended (21g) threshold replaced by 0.085
App stopped. Returning the thresholds & cells selected.

runSCENIC_4_aucell_binarize(scenicOptions)

Binary regulon activity: 172 TF regulons x 9537 cells.
(172 regulons including 'extended' versions)
154 regulons are active in more than 1% (95.37) cells.

ds1

fileloc <- scenicOptions@fileNames

temp <- lapply(fileloc[["output"]],lamb,"ds1_SCENIC/") %>% as.character() %>% as.matrix()
rownames(temp) <- rownames(fileloc[["output"]])
colnames(temp) <- "fileName"
fileloc[["output"]] <- temp

temp <- lapply(fileloc[["int"]],lamb,"ds1_SCENIC/") %>% as.character() %>% as.matrix()
rownames(temp) <- rownames(fileloc[["int"]])
colnames(temp) <- "fileName"
fileloc[["int"]] <- temp

scenicOptions@fileNames <- fileloc
scenicOptions@settings[["tSNE_filePrefix"]] <- "ds1_SCENIC/int/tSNE"

regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
regulonActivity_byCellType <- sapply(
  split(rownames(cellInfo), cellInfo$CellType),
  function(cells) rowMeans(getAUC(regulonAUC)[, cells])
)

regulonActivity_byCellType_Scaled <- t(scale(t(regulonActivity_byCellType), center = T, scale = T))
topRegulators <- reshape2::melt(regulonActivity_byCellType_Scaled)
colnames(topRegulators) <- c("Regulon", "seuratCluster", "RelativeActivity")
topRegulators <- topRegulators[which(topRegulators$RelativeActivity > 0), ]

aucellApp <- plotTsne_AUCellApp(scenicOptions, logMat)

savedSelections <- shiny::runApp(aucellApp)


plotTsne_AUCellHtml(scenicOptions, logMat,"ds1_scenic")

newThresholds <- savedSelections$thresholds
scenicOptions@fileNames$int["aucell_thresholds", 1] <- "ds1_SCENIC/int/newThresholds2.Rds"
saveRDS(newThresholds, file = getIntName(scenicOptions, "aucell_thresholds"))
saveRDS(scenicOptions, file = "ds1_SCENIC/int/scenicOptions.Rds")
#save.image()

runSCENIC_4_aucell_binarize(scenicOptions)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.